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We present three dimensional realizations of the model introduced recently by (Karlin, Bosch, 
Chikatamarla, Phys. Rev. E 2014; Ref m) and review the role of the entropic stabilizer. The 
presented models achieve outstanding numerical stability in presence of turbulent high Reynolds 
number flows. We report accurate results for low order moments for homogeneous isotropic decaying 
turbulence and second order grid convergence for most assessed statistical quantities. The explicit 
and efficient nature of the scheme renders it a very promising candidate for both engineering and 
scientific purposes in the vicinity of highly turbulent flows. 
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I. INTRODUCTION 

The lattice Boltzmann method (LB) [21 [3] is a mod¬ 
ern and highly successful kinetic-theory approach to 
computational fluid dynamics (CFD) and computational 
physics of complex flows and fluids, with applications 
ranging from turbulence |4] to flows at a micron scale |5] 
and multiphase flows iniizi, to relativistic hydrodynamics 
[5], soft-glassy systems jS] and beyond. 

While conventional methods for CFD attempt to solve 
the Navier-Stokes equations, LB’s underlying equations 
form a kinetic system which can be thought of as a dis¬ 
crete analogue to Boltzmann’s equation. Carefully de¬ 
signed populations corresponding to a set of dis¬ 

crete velocity vectors Vi, i = 1,... ,b spanning a regularly 
spaced lattice with nodes x. Conceptually, the dynamics 
of populations fi can be split into free flight (advection) 
and collison (relaxation) which is reflected by the general 
one-parametric LB equation 

f,{x + v„t+l) = fi = (l-/3)/i(a;,t)-b/?/“'"(a;,t). (1) 

Here the left-hand side is the propagation of the popu¬ 
lations along the lattice links, while the right-hand side 
is the so-called post-collision state /'. The relaxation 
parameter fj is associated with the transport coefficient 
of the macroscopic target equation (the incompressible 
Navier-Stokes equations herein). The mirror state 
represents the maximally over-relaxed state. Realization 
of hydrodynamics was made possible, in the hrst place, 
with the lattice Bhatnagar-Gross-Krook (LBGK) model 
[ini[n], in which one takes 

/r" = 2fp - /,. (2) 

Here is the equilibrium which is found as a maximizer 
of the entropy El. 

S[/l = -i:/.ln(|-). (3) 
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subject to hxed locally conserved helds identified by the 
first D + 1 velocity moments, p = fi (density) 

and pu = X]i=i(momentum density), and where 
the weights Wi are lattice-specific constants. With the 
proper symmetry of the lattice, the LBGK equation, 0 
and ([^, recovers the Navier-Stokes equation for the fluid 
velocity u, with the kinematic viscosity 


where Cg is the speed of sound (a lattice dependent 
0(1) constant). The LBGK model is unambiguous since 
P € (0,1) is fixed by the kinematic viscosity Q. The 
limit /I —)■ 1 (small kinematic viscosity) is particularly 
important as it is pertinent to achieving, if only in princi¬ 
ple, high Reynolds number regimes. A notable departure 
from the continuous BGK approximation of Boltzmann’s 
equation becomes manifest in the feature of over-relaxing 
towards the mirror-state which disconnects LBGK from 
the kinetic theory domain /3 G (0, [IB] . 

Almost immediately after its inception, the LBGK 
model has taken lead in the lattice Boltzmann approach 
to the simulation of complex hydrodynamic phenomena 
mm, and remains the “working horse” of the LB meth¬ 
ods to-date. Popularity of LBGK is primarily based on 
its simplicity and exceptional computational efficiency. 
Despite of its promising nature and popularity the LBGK 
model shows severe dehciencies (disruptive numerical in¬ 
stabilities) already at moderate Reynolds numbers. This 
precluded the LB method from making a sustainable im¬ 
pact in the field of computational fluid dynamics. 

A number of approaches can be found in the literature 
intended to alleviate this issue. We will restrict the fol¬ 
lowing short discussion to methods without explicit tur¬ 
bulence models. Most notably, the entropic lattice Boltz¬ 
mann method (ELBM) features non-linear stability and 
has shown excellent performance [T5lll7j . While ELBM 
converges to LBGK in the resolved case, it locally alters 
the relaxation parameter which in turn modihes the vis¬ 
cosity in order to fulfil the second law of thermodynamics 
by both enhancing and smoothing the features of the flow 
where necessary, subject to an entropy condition. 
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The dimension b of the kinetic space populated by fi 
is usually greater than strictly necessary for recovering 
the Navier-Stokes equations. In three dimensions, for 
instance, ten linearly independent basis vectors would 
suffice to represent the conserved variables and the sym¬ 
metric stress tensor liap = Although it 

can be shown that the coupling to the non-hydrodynamic 
higher (’’ghost”) moments cannot be chosen arbitrarily 
in the limit Ma —>■ 0 [TB] independent relaxation of these 
moments may have a favourable effect on the numerical 
stability. 

The multiple relaxation parameters (MRT) LB 
schemes follow this line of thinking. While the relaxation 
of the off-diagonal parts of the stress tensor are fixed 
by the choice of kinematic viscosity, the MRT scheme 
embraces that the relaxation of higher order moments 
should not affect the dynamics of the flow field (up to the 
Navier-Stokes level) and hence can be used to construct 
more stable LB schemes. Based on a separation of scales 
between fast and slowly varying moments several MRT 
schemes were suggested for the choice of relaxation of 
higher order moments (beyond the stress tensor) [TBIEB] . 
The choice of relaxation parameters is crucial in order to 
increase the operational range in terms of stability and 
requires careful tuning. 

In order to eliminate the influence of higher order mo¬ 
ments which may be rapidly oscillating and cause numeri¬ 
cal instabilities the regularized LB scheme (RLB) 
was proposed. While the relaxation parameters for MRT 
schemes are tunable, they are fixed for RLB such that 
the higher order non-equilibrium moments vanish. As 
any finite lattice representation introduces discrete arti¬ 
facts among the higher order moment tensors, the regu¬ 
larization operation ensures isotropy, albeit in the con¬ 
fined subspace limited up to the stress tensor level. Al¬ 
though MRT and RLB models were successful in slightly 
stabilizing the LB method, they still remain challenged 
by high Reynolds numbers [25]. 

Recently, the authors have developed a scheme with¬ 
out the need for tunable parameters or turbulent viscos¬ 
ity (Karlin, Bosch, Chikatamarla, Phys. Rev. E 2014; 
Ref HI) which has demonstrated a significant extension 
in the operation range for simulations at high Reynolds 
numbers. Promising results have been reported for both 
two and three dimensions, as well as for complex bound¬ 
aries and in presence of turbulence. Much like ELBM, 
entropic considerations have been employed to render the 
scheme stable without introducing considerable compu¬ 
tational overhead and by keeping the simplicity and lo¬ 
cality of the LBGK and MRT schemes. Below we shall 
refer to this class of models as KBC models for brevity. 

In the remainder of this manuscript we present the 
construction of a family of KBC models for the standard 
D3Q27 lattice and review the details of the entropic sta¬ 
bilization. The Kida vortex flow and a randomly gen¬ 
erated initial condition shall be employed as turbulent 
benchmark flows and are studied in detail for one of the 
realizations of KBC and the stability domain of LBCK, 


RLB and KBC variations is assessed numerically. 


II. EQUILIBRIUM ON THE STANDARD 
LATTICE IN THREE DIMENSIONS 

In this section, we shall review the standard lattice 
in three dimensions and the corresponding equilibrium. 
While the material presented in this section is of a review 
character, we shall highlight some fundamental features 
of the equilibrium which were not fully discussed so far 
in the literature. 

The standard lattice in a dimension D is built as ten¬ 
sor (direct) product of D copies of the fundamental one¬ 
dimensional set = i, where j = 0, ±1: 

v-i = -1, vq = 0, vi = 1. (5) 

The natural il)-dimensional Cartesian reference frame 
generated by the tensor product of D copies of the fun¬ 
damental set ([^ makes it convenient to enumerate the 
discrete velocities accordingly. Considering the three- 
dimensional case D = 3 below, we write for any of the 
b = 27 discrete velocities, 

Vi — (Vix^ ‘^iy, ^ — I 5 ■ • • ? 27, Via ^ (5) 

The equilibrium on this standard product-lattice max¬ 
imizes the entropy (|^ subject to fixed conservation laws 
of density and momentum m- It is written most ele¬ 
gantly in the following product-form: 


[Biux)r- [BiuyT^ [BMT- , 


where the weights Wi are 


(7) 


W, = (8) 

and function 'k is universal for all the discrete velocities 
(it does not depend on the discrete velocity index). 


^'(m) = A{ux)A{uy)A{u^), (9) 

with 

A{u) = 2- ^/l + 3u^. (10) 


Furthermore, the universal function B{u) contributing to 
the formation of the last term in Q is written as 


B{u) 


2u + Vl + 
1 — u 


( 11 ) 


Finally, the weights Wa in Q are dictated by the speed 
of sound, 


Cs = 


1 

71 ’ 


( 12 ) 


and are 


Wo = 2/3, W-i = 1/6, Wi = 1/6. 


(13) 
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Remind that the Maxwellian in D dimensions is writ¬ 
ten as a product with respect to an arbitrarily fixed 
Cartesian reference frame, in accord with the familiar 
property of the shifted Gaussian distribution, 



(14) 


It is easy to recognize the Maxwellian character of the 
product-lattice equilibrium. The multiplication of the 
weights in Q corresponds to the first multiplier, the 
function 4' corresponds to the second multiplier, while 
the product of functions B reflects the last multiplier. 

However, the true Maxwellian is isotropic, as also re¬ 
flected by reading equation ( [T4| ) from right to left, the 
products collapse to a dependence on the kinetic energy 
in the co-moving frame alone, and the reference to the 
arbitrarily fixed Cartesian coordinates disappears. This 
is not so with the discrete velocities. It is imperative 
therefore to demonstrate that the product-form Q is 
manifestly isotropic to the order of accuracy of the lat¬ 
tice Boltzmann model. This can be done most elegantly 
in the following way: Instead of expanding each popula¬ 
tion Q into powers of the velocity components Ua , let us 
first expand the logarithm of (we consider a generic 
case of D below): 


truncation of the expansion of any function ip around 
M = 0 to get 


[ln/r ]2 = Inp+lnlTi-f ^1 - ^(m • u)^ +3{v,-u), (16) 
where we have used the standard notation, 

D 

CL • h — ^ ^ 

Q — 1 

for the Cartesian scalar product of D-dimensional vec¬ 
tors. Then, using the identity, [fP ]2 = [exp([ln/®‘^] 2 )] 2 , 
we get 


{f:\ = 


U ■ Vi 


{u ■ Vi)^ — Cs(m • u) 


(17) 

The second-order polynomial ( [T7| ) generated by the equi¬ 
librium 0 is manifestly isotropic, a nd w ith the definition 
of the speed of sound Cg = l/-\/3 (12) it is identical to 
the standard lattice Boltzmann equilibrium. 

III. MOMENT SYSTEM FOR D3Q27 


We consider the standard 27 velocity lattice (D3Q27). 
Let the natural moments be defined as 


D D 

InfP = lnp + InITi + ^ InB(ua) + ^ Via \nB{ua). 

Q!=l 

(15) 

Let us denote [p{u )\2 the operation of the second-order 


pMpqr = {fivtvIyvQ , p,q,r e {0,1, 2} (18) 

where notation (...) is used as a shorthand for summation 
over all velocity indices. In the basis spanned by the 
natural moments populations can be represented as 


^<^Qxzz + 3crAfi22 ~ 3117202 ~ 31172 2 0 + 3117222 ) 
BAQyzz + 3 XM 212 — 3117o22 ~ 3117220 + 3117222) 

'Qyyz + 3(5117221 — 3Mq22 — 3117202 + 3117222) 
A 7 i 22 ~ AlV/212 ~ (tA 117 ii 2 — M222) 

B122 — SM221 — cr 5 Mi 2 i — M222) 

'I212 ~ SM221 — XSM211 — M222) 
b (tA 117 ii 2 -f a'6Mi2i X6M211 -b M222) ■ 


(19) 


f{0,0,0) — p{^ ~ T + A7o22 + M202 + M220 — M222) 
f{(7,0,0) ~ qP (3(TUa, -b ‘INxz ^yz “b T ^(jQxyy 
f{o,x,o) = ^P (3Aity — Nxz + ^Nyz + T — 3XQxxy — 

f{0,0,S) = \p (SSUz: — Nxz — Nyz: + T — iSQxxz 
f{<7,\,0) — \p{<xXIVxy -b XQxxy + CrQxyy + M22O 
f{a,0,5) = ip(o-(5n xz + dQ 

xxz + ctQ xzz + 717202 
f{0,X,S) 4P {XSTly^ -b SQyyz XQyzZ T X(Iq 22 

f{a,X,6) = ^P {(xXSQxyz + 0 'Mi 22 + XM212 + SM22I 

Here we depart form the usual single index i and use a 
somewhat unconventional indexing for the discrete veloc¬ 
ities, where indices cr, A ,7 G {—1,1}. Note that labelling 
of the velocities by a triad (•,•,•) is unambiguous as long 
as the first, the second, and the third entries are always 
associated with the x, the y, and the z coordinates, re¬ 
spectively, in the once fixed Cartesian frame. 


We chose to rename some of the natural moments as 
reminiscence to their physical meaning: 

T = M 200 + M 020 + M 002 

is the trace of the stress tensor at unit density, 

Nxz = M200 — M002 
Nyz = Mq2Q — Mqq2 
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are the normal stress differences at unit density and 


^xy — -^110 
Ilxz = Mioi 
Ilyz = Mon 

are the off-diagonal components of the stress tensor at 
unit density. The third order moments lack a direct 
physical interpretation in the isothermal case but are de¬ 
noted as Qxzz — .M1025 Qxxy — .M21O; Qyyz — M021, 
Qxxz — Tf 20 l 5 Qyyz — M021 and Qxyz — ^ 111 - 

Another basis is given by the central moments of the 
form 


P^pqr — ifiip^ix ‘^x)^ {'^iy {piz ^ 2 ) ) • (^ 0 ) 


Using the mapping from natural to central moments 
(A.l)-(A.23), which is linear in the non-conserved mo¬ 
ments, a similar moment representation of the popula¬ 
tions in the central moments basis can then be writ¬ 
ten (|A.24|)-(|A 311). 


IV. KBC 

Let us then group the contribution to each population 
in three parts 


fi — ki -\- Si hi^ 


( 21 ) 


where ki (= kinematic part) depends only on the locally 
conserved fields, Si (= shear part) depends on the (devia- 
toric) stress tensor Ho and hi (= higher-order moments) 
is a linear combination of the remaining higher-order mo¬ 
ments. Representation (21 1 is easily obtained for any 


lattice and any moment basis. 

With the representation (21 1 , a different mirror state 


can be sought in a one-parameter form, 

/r" = h + [25^ - Si] + [(1 - 7)/i* + ihTl (22) 

where 7 is a parameter which is not yet specified and the 
terms s®*^ and denote the s and h part evaluated at 
equilibrium 0- When ( [ 2 ^ is used in 0 , one arrives 
at nothing but a special (not the most general) MRT 
model. For any 7 , the resulting LB model still recovers 
hydrodynamics with the same kinematic viscosity v 0 . 
Note that Eq. 0 and Eq. (p^ can be also rewritten as 


/' = /, + 2/3(/^J^ - U) 


with the generalized equilibrium [26l - [^ of the form, 

/f'' = /r+(i/2)(7-2)(hr-/io. 


For 7 = 2 we obtain the LBGK model while the choice 
of 7 = ^ results in the generalized family of regularized 
LB (RLB) models. 

While the dependence of s on the deviatoric stress ten¬ 
sor no is mandated, one may include further moments. 


index k/f. 

’ d/p 


t/p 


<iIp 


(0,0,0) 1 

0 


-T 


0 


(cr, 0,0) ctm 3;/2 {2N^z-Nyz 

)/6 

T/6 

^(^Qxyy H” Qxz: 

0/2 

(0,A,0) Xuyf2 (—+ 2A^y, 

.)/6 

r/6 

xxy “1” Qyz 2 

0/2 

( 0 , 0 , (5) 5uz/2 {—Nxz — Nyz 

)/6 

r/6 

—S{Qxz 

“h Qyyz 

0/2 

(a,A,0) 0 

aXIlj:y/4: 


0 

xxy 

+ O’Qxyi 

0/4 

(a, 0,5) 0 

a&kdzz/^ 


0 

(^^Qxxz 

O'Qxzz 

)/4 

(0,A,5) 0 

A5n„./4 


0 

{SQyyz 

“h XQyzz 

)/4 

(ct, A,5) 0 

0 


0 

(JX^^^xyz /8 



Table I. Moment groups in the natural basis listed for all 
discrete velocity directions, k is the kinematic part, d is a 
function of the moments amounting to the deviatoric stress 
tensor, t is a function of T = tr(II) associated with the fluid 
bulk viscosity and q contains the contributions of the third 
order moment rank-3 tensor. 


provided the basic symmetry properties are not violated. 
We will consider eight different KBC models, all of which 
give the correct kinematic viscosity in the hydrodynamic 
limit but differ in the choice of Si. 

The four natural moment KBC models considered 
herein are characterized by 

s G {d, d-t-t, d-t-Q', d t , 


where d depends on the deviatoric stress tensor Ho, t 
depends on the trace of the stress tensor T = tr(n) and 
q is a function of the third order moments Qap-y- Table ^ 
shows a comprehensive list for d, t, q and the kinematic 
part k in the natural moments representation. Note that 
the model with s = d + t and fixed 7 = ^ is equivalent 
to the regularized LB model in [21]. 

Analogously, we define four KBC variations in the cen¬ 
tral moment basis, 

s G \ d^ d-t-t, d q^ d t 


where fc, d, t and q correspond to the central moments 
representations (A.24)-( A.31). The remaining higher or¬ 
der parts are trivially given as 


h = f — k — s 


in all cases. 

The major change of perspective here is that the sta¬ 
bilizer 7 should not be considered as a tunable param¬ 
eter. Rather, it has to be put under entropy control 
and computed by maximizing the entropy in the post¬ 
collision state /'. This matches the physics of the prob¬ 
lem at hand, since constrained equilibria correspond to 
the maximum of the entropy (here the constraint is 
that the s part remains fixed by the over-relaxation, 

gmirr ^ 2sf - s,). 

Specifically, let «S'( 7 ) be the entropy of the post¬ 
collision states appearing in the right hand side of 0 , 
with the mirror state ( |22| . Then we require that the sta¬ 
bilizer 7 corresponds to maximum of this function. In¬ 
troducing deviations, Asi = Si — s'p and Ahi = hi — h®'^. 
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the condition for the critical point reads: 


b 

E 

i=l 


Ahi In ( 1 + 


(1 - /3^)Ahi - (2/3 - l)Asi 


= 0 . 


(23) 

Equation (23) suggests that among all non-equilibrium 
states with the fixed mirror values s™“'' = 2 s®'^ — Si, we 
pick the one which maximizes the entropy. In contrast 
to MRT, the entropic stabilizer 7 is not tunable but is 
rather computed at each lattice site in every time step 
from Eq. (23). Thus, the entropic stabilizer self-adapts 


to a value given by the maximum entropy condition (231. 


In order to clarify the properties of the solution to Eq. 


(23), let us introduce the entropic scalar product 
in the ^-dimensional vector space. 


(w = E^. 

Ji 


(24) 


and expand in (23) to the first non-vanishing order in 


V. NUMERICAL SIMULATIONS 

The Kida vortex flow is a well-studied benchmark flow 
which evolves from a simple deterministic and symmet¬ 
ric initial condition to a state which resembles a fully 
developed turbulent flow, which features a correspond¬ 
ing energy cascade. The initial conditions for the flow 
field are given by 

Ux{x^ y, z) = Uq sina;(cos 3j/cos z — cosy cos 3^;) 

Uy{x, y, z) = Uq siny(cos 3z cosx — cos 2 : cos 3a;) (26) 

Uz(x, y, z) = Uq sin 2 :(cos 3 a; cosy — cos a; cos 3y) 

where x,y,z G [ 0 , 27r] and periodic boundary conditions 
are imposed in all directions. The Reynolds number is 
defined as Re = UqN/ v where N is the domain size. Ini¬ 
tial conditions for the density (and pressure p = pc^) 
and higher order moments are obtained by solving the 
convection-diffusion equation + V • (pwo) = D Ap on 
the same grid beforehand until steady state is reached in 
a similar process as described by [53]. 

The Kida vortex flow has been analyzed extensively 
using DNS [501155] . The evolution of enstrophy shows a 
steep increase in the early stage of the simulation and 
reaches a maximum value before it decays. For the con¬ 
vergence study we investigate data collected from time 
points around the peak of enstrophy which indicates the 


Asi /and Ahi /to obtain 


!_/ 1\ (A^ 

^ \ p) {Ah\Ahy 


(25) 


The result (25) explains the mechanism of failure of the 
proposal 7 « 1 at /3 « 1: Whenever vectors As and Ah 
are non-orthogonal (in the sense of the entropic scalar 
product), the deviation of 7 * from 7 = 1 may become 


very significant. Indeed, in (25), the correlation between 
the shear and the higher-order parts ~ {As\Ah) is not a 
correction to 7 = 1 but rather a contribution of same 
order 0 ( 1 ). 


We have found that the estimate (25) was sufficient for 


stabilizing all the simulations which renders the KBC an 
explicit and efficient method with only slightly increased 
computational costs compared to the standard LBGK. 
The resulting collision operation is given here: 

1 compute conserved quantities p^Ua 

2 evaluate equilibrium f/'^{p,Ua) 

3 compute s and s®'^, i.e. s = d + t + q, 

-t (see table 

4 compute Asi = Si — 

5 compute Ahi = hi — h^ = fi — jp — Asi 

6 evaluate 7 from equation (25) 

7 relax f/= fi - P {2Asi + ^Ahi) 


existence of large gradients which are often numerically 
challenging. A simulation was considered stable if it run 
until the mean enstrophy 

was sufficiently decayed (D/Dq < 5%) and 
Ua =Ua- (Ua) . 

In order to assess the stability region, the domain size 
N = 100 and initial velocity Uq = 0.05 were fixed and 
the Reynolds number Re was increased in steps of 500. 
While LBGK seized to yield sensible values at Re > 5000, 
ELBM was always stable (tested up to Re = 10^). Like¬ 
wise, all KBG models were always stable, independently 
of the moment basis or the choice of s. The only RLB 
models, however, which achieved the same stability are 
using s = d and s = d, respectively. All other six RLB 
models, among them the standard model s = d + t, were 
less stable than LBGK. 

Accuracy of KBC scheme is studied in detail using 
Kida vortex flow at Re = 6000. KBC model s = d + t + q 
is used exemplarily in all further simulations and com¬ 
pared to LBGK simulation at N = 600 (run D) where the 
flow is considered to be resolved as indicated by the Kol¬ 
mogorov length scale p = {p ~ 1.2 lattice units 
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t = 0.25 



t = 0.5 




t 

= 0.75 



A 

B 

C 

D 

P 

A 

B 

C 

D 

P 

A 

B 

C 

D 

P 

k- 10^ 

8.528 

8.644 

8.663 

8.657 

2.21 

6.237 

6.385 

6.411 

6.402 

2.10 

3.808 

3.947 

3.965 

3.954 

1.87 

OA^ 

1.472 

1.735 

1.784 

1.784 

2.67 

1.919 

2.786 

3.173 

3.216 

2.46 

1.551 

2.130 

2.262 

2.301 

2.13 

eA • 10® 

2.459 

2.892 

2.974 

2.974 

2.65 

3.247 

4.646 

5.288 

5.360 

2.44 

2.646 

3.553 

3.770 

3.835 

2.10 

S 3 ■ 10" 

10.73 

7.361 

9.071 

10.34 

1.23* 

28.91 

38.86 

47.78 

48.97 

2.04 

29.37 

33.81 

39.71 

39.99 

2.62 

54 

5.075 

6.188 

6.176 

6.149 

2.66 

3.993 

5.218 

6.063 

6.159 

2.25 

4.140 

4.319 

4.840 

4.842 

4.23 

55 

2.902 

2.157 

2.472 

2.891 

0.81* 

3.210 

6.084 

9.044 

9.690 

1.66 

4.337 

4.313 

6.074 

5.978 

2.05 

56 

52.50 

83.87 

81.18 

79.34 

1.93 

34.07 

66.86 

96.83 

100.1 

2.17 

38.64 

36.27 

51.73 

51.55 

3.08 

lo/N ■ 10' 

10.13 

8.788 

8.573 

8.564 

3.72 

4.796 

3.473 

3.070 

3.022 

2.60 

2.809 

2.207 

2.094 

2.050 

2.05 

uo ■ 10^ 

2.920 

2.940 

2.943 

2.942 

2.23 

2.497 

2.527 

2.532 

2.530 

2.02 

1.951 

1.987 

1.991 

1.988 

1.81 

ro/N 

34.68 

29.89 

29.13 

29.11 

4.06 

19.21 

13.74 

12.12 

11.94 

2.67 

14.39 

11.11 

10.52 

10.31 

2.14 

Reo 

3549 

3101 

3028 

3024 

3.52 

1437 

1053 

932.7 

917.6 

2.55 

657.8 

526.2 

500.3 

489.2 

1.96 

A/A • 10^ 

5.376 

4.991 

4.927 

4.925 

3.91 

4.001 

3.384 

3.179 

3.155 

2.57 

3.463 

3.043 

2.960 

2.931 

2.10 

ux ■ 10^ 

2.384 

2.401 

2.403 

2.402 

2.08 

2.039 

2.063 

2.067 

2.066 

2.38 

1.593 

1.622 

1.626 

1.624 

1.98 

rx/N 

2.254 

2.079 

2.050 

2.050 

2.81 

1.962 

1.640 

1.537 

1.527 

2.72 

2.174 

1.876 

1.821 

1.805 

2.26 

Rex 

153.8 

143.8 

142.1 

142.0 

3.44 

97.89 

83.79 

78.86 

78.21 

2.46 

66.22 

59.23 

57.75 

57.11 

1.92 

p/N ■ 10® 

2.202 

2.115 

2.100 

2.100 

2.77 

2.055 

1.879 

1.819 

1.813 

2.67 

2.163 

2.009 

1.979 

1.971 

2.29 

Ur, ■ 10® 

3.784 

3.940 

3.968 

3.968 

2.72 

4.056 

4.436 

4.582 

4.597 

2.59 

3.853 

4.148 

4.210 

4.228 

2.19 

Tn/N ■ 10® 

5.821 

5.368 

5.293 

5.293 

2.82 

5.066 

4.235 

3.970 

3.943 

2.69 

5.612 

4.843 

4.701 

4.662 

2.30 


Table II. Comparison of LBGK and KBC {s = d + t + q) for statistical quantities in Kida vortex flow at Re = = 6000 

and t = ^ — 0.25, 0.5, 0.75. Resolutions N = 100, 200,400,600 for KBC runs A, B, C and resolved LBGK run D, respectively. 
Convergence rate p of error w.r.t. LBGK solution estimated from polynomial fit (* indicates exclusion of lowest resolution). 
All gradient-based quantities are computed with eighth-order of accuracy central difference stencils. 

Turbulence characteristics: length scale Zq = velocity scale Uq = time scale tq = Iq/uq, and Reynolds number 

Reo = Zowo/i^. 

Taylor characteristics: Taylor micro scale A = velocity scale u\ = u' = (2fc/3)^'^^ (r.m.s. turbulence intensity), 

time scale tx = X/ux and Reynolds number Re^ = Xuxjv- 

Kolmogorov characteristics: length scale p = velocity scale Un = and time scale r,, = 



1.25e-02 2.50e-02 3.75e-02 5.00e-02 6.25e-02 


Figure 1. Iso-surface of vorticity component = 0 at time 
t = 0.5 colored with velocity magnitude rendered at 2 ; = 0, 
x,y € [0,-7r] plane. Runs A {N — 100), B (A = 200), C 
{N = 400) and D (A = 600, reference solution). 


where 

^ + Sif) + §71)) 

is the rate of turbulence kinetic energy dissipation. Reso¬ 
lutions N = 100, 200,400 are considered in the following 
(runs A, B and C, respectively). Convergence towards 
resolved LBGK simulation is reported in table |ll] for var¬ 
ious statistical quantities. Unless stated otherwise all 
quantities are given in lattice units. 

Figure shows a comparison of the vortex structures 
for the four simulations roughly at the point of maximum 
enstrophy. The vorticity configuration is thus affected by 
the under-representation of large gradients in the coarse 
resolntions (see next section). Nevertheless, the large 
vortex structures are akin at all resolutions. The largest 
KBC simulation (run C) is hardly discriminable form the 
reference LBGK simulaion (run D). 

A. One-Point Statistics 

Enstrophy O and turbulence kinetic energy 
k = ^ («) 

are important global quantities characterizing the flow 
and its history. Figures and respectively, depict the 
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Figure 2. Mean enstrophy evolution with time for simulations 

A (— ■ ■), B ( ■), C (-) and D (—). Gradients evaluated 

with second-order of accuracy. 



Figure 3. Kinetic energy evolution with time for simulations 
A(- - ). B( ■), C(—)andD(-). 


evolution of both quantities with time t = N/Uq. It 
is apparent that in the under-resolved KBC simulations 
N = 100, 200 the enstrophy peak values are not well rep¬ 
resented. However, for coarse resolutions this is expected. 
The kinetic energy on the other hand decays quite sim¬ 
ilarly for all simulations. Table reports the numbers 
at three selected time instances. During simulation, gra¬ 
dients are evaluated using second-order finite differences, 
which are solely used for reporting the enstrophy evo¬ 
lution in figure All quantities based on gradients in 
table [nj however, are computed with a nine point cen¬ 
tral difference stencil of eighth-order accuracy, which also 
explains the discrepancies in enstrophy between figure 
and table El for the lowest resolution N = 100. 

While the energy seems to be dissipated similarly with 
time it is important to study the kinetic energy k and dis¬ 
sipation rate e thereof across flow scales in order to decide 
whether low-order statistics of turbulent flows yield sen¬ 
sible values in coarse resolution simulations despite the 
under-representation of high gradients. The instrument 
at hand is the spectral representation of the kinetic en¬ 
ergy distribution E{k) where k = ||k|| is the modulus 


10 ^ 

10 " 

10-2 
10 -‘ 

10 -" 

^ 10 -" 

2 10 - 1 " 

&q 10-12 
10 -1* 

10-16 
10 - 1 ® 

10 - 2 " , „ , 
10-1 10" lOi 

KT] 

Figure 4. Kinetic energy density spectrum at t = 0.5 for 

simulations A (— ■ ■), B ( ■), C (-) and D (—). The dotted 

line with slope —5/3 shows the Kolmogorov scaling. 



0.5 1.0 1.5 2.0 


Figure 5. Cumulative distribution function of the energy- 
dissipation rate density D{k) at time t = 0.5 for simulations 
A(- - ). B( ■), C(—)andD(-). 

of the wave number vector. Figure shows the non- 
dimensional energy density distribution normalized with 
kinetic energy 

pOO 

k= E{k) (Ik. 

Jo 

According to [331 ES] the energy scales as E in 

the inertial sub-range. The studied Kida flow here does 
not exhibit large enough Reynolds numbers to see an ex¬ 
tended inertial range. However, it is apparent that the 
energy scales similarly across resolutions and a sharp cut¬ 
off is visible at the smallest scales. This indicates that the 
KBC model is capable of producing the expected energy 
distribution throughout the scales without an explicit 
turbulence model. A case for higher Reynolds number 
shall be examined below. 

The cumulative distribution function of the energy- 
dissipation rate density 

D{k) = 2vkJE{k) 





















Figure 6 . Second order longitudinal structure function at t = 

0.5 for simulations A (— ■ ■), B (—■), C (-) and D (—). The 

dotted line indicates the theoretical scaling. 



Figure 7. Longitudinal and transversal velocity correlation 

functions at t = 0.5 for simulations A (— ■ ■), B ( ■), C (-) 

and D (—). 


where 


e = / D(k) dK 

Jo 

illustrates the scales of eddies responsible for the dissipa¬ 
tive process, see figure The under-resolved simulations 
employ expectedly larger eddies for the bulk of the dissi¬ 
pation (see also table [iT] where e is reported). 

The longitudinal skewness factor 



is another global statistical quantity in real space which 
we report in table |llj In agreement with figure we find 
that the outcome of the lowest resolution fV = 100 is 
rather inconsistent with the trend observed in the other 
simulations, however, agrees well with the resolved case. 
The lower convergence rate for the odd-order skewness 
factors may be caused by the inherent lack of isotropy in 
the third-order moments. However, further studies are 
needed to draw a concise conclusion. 

The remainder of table [H] is a compilation of the tur¬ 
bulence, Taylor and Kolmogorov flow scales. Here and 
with the vast majority of the reported quantities we ob¬ 
serve a grid convergence rate of « 2 as expected in the 
context of LB simulations. 


B. Two-Point Statistics 


The longitudinal structure function of order n defined 
as 

exhibit linear scaling on logarithmic plots [Ml[55] . In par¬ 
ticular, the second order structure function scales with 


Bfi ^ Figure [^depicts the results with the theoret¬ 

ical scaling. Due to the relatively low Reynolds number 
we may not identify an extended inertial range but we 
note that the simulations agree well with the reference 
over the entire range of r. 

Another real-space two-point statistical quantity that 
can be used to assess different numerical techniques is the 
correlation of the velocity field. Here the longitudinal and 
transversal correlation functions are defined as 


n / ^ _ {Kix,y,zy^{x + r,y,z)) 

r11V / III \ ! { \\ 

{K{x,y,zy{x,y,z)) 
n , ^ y, z)u' {x + T, y, z)) 

pUr) = -frT- 

{u'y{x,y,z)u' {x,y,z)} 


A comparison at time t = 0.5 is given in figure All 
simulations with N > 200 show excellent agreement with 
the reference solution. At the maximum distance r = 
0.5 the velocity components are still correlated which is 


associated with the low Reynolds number (see figure 10 
for comparison). 



Figure 8 . Evolution of the mean entropic stabilizer { 7 ) and 
its standard deviation with time for simulation A. 
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Figure 9. Turbulence kinetic energy spectrum at Re^ = 580, 
t = 0.75 (solid) and initial spectrum at Rca = 2000 (dashed). 
The dotted line indicates the theoretical scaling. 



Figure 10. Longitudinal and transversal velocity correlation 
functions at Rca = 580, t — 0.75 (solid) and initial values at 
Rex = 2000 (dashed). 


C. Stabilizer 7 

The importance of the self-adjusting stabilizer 7 be¬ 
comes clear when considering the evolution thereof with 
time. Figure]^ shows the mean stabilizer ( 7 ) at TV = 100 
with its standard deviation. It is apparent that 7 is far 
form being constant. In fact, its evolution is closely cor¬ 
related with the flow state. In regions of high turbulence 
intensity it is distinctly different in the mean and shows 
larger fluctuations. In [1] we also report the distribution 
of 7 in space which again reveals the close relation of 7 
to the flow. It is therefore conjectured, and supported by 
our simulations and comparisons to literature, that any 
MRT using 7 = const may not reach numerical stability 
to the same extent. 


D. Large Reynolds Numbers 

The second numerical example considered in this paper 
is the simulation of a highly turbulent flow starting from 


a random initial condition and decaying with time. The 
periodic and cubic domain with box-length N = 400 was 
initialized with a flow field generated from a prescribed 
narrow-banded initial energy spectrum peaked at grid 
wave-number « 8 , 

i?o = 400 ^ 2 ^Uo^/ 7 r^ exp 

5 = 20000 (f)^^%/wo, 

where n = 8.164970 • 10“® and Uq = 0.01. The initial 
values of the density (and pressure) and the higher order 
moments were generated with the identical procedure as 
described earlier for the Kida flow. 

The main objective is to test the KBC scheme for an 
under-resolved simulation (Kolmogorov length 77 « 10 
lattice units) at large Reynolds numbers (Re^ « 600). 
In particular, we ask whether the scheme is stable for a 
random and highly turbulent flow in absence of a deter¬ 
ministic and highly symmetric initial condition, whether 
low-order statistics are well represented and physical dis¬ 
sipation (i.e. scaling laws) is modelled correctly. By 
means of this simulation we examine the general ques¬ 
tion of the performance for large Reynolds numbers in 
an under-resolved simulation from yet another point of 
view. While a resolved simulation was not attempted, we 
compare our results to the classical scaling laws. 

Figurej^shows the turbulence kinetic energy spectrum 
at t = N/uq = 0.75. The inertial range is extended 
and the scaling is more apparent than in the less turbu¬ 
lent simulations above. The sharp cut-off at the smallest 
scales is still maintained despite the coarse resolution. 
As before, numerical stability is naturally guaranteed to 
very high Reynolds numbers. While the initial spectrum 
is narrow and steep, it flattens during the course of en¬ 
ergy decay and exhibits the Kolmogorov scaling in the 
inertial sub-range roughly at peak of mean enstrophy. 
Figure [l0| shows the velocity correlations where the con¬ 
tributions to the correlations are vanishing for r/N > 0.2 
at t = 0.75. Hence, the velocity field is largely uncorre¬ 
lated as one would expect from isotropic homogeneous 
turbulent flows. While these results are far from a com¬ 
prehensive study, they contribute to the overall assess¬ 
ment that the KBC scheme might perform well even in 
the case of severe under-resolution. A more comprehen¬ 
sive investigation shall be conducted in a further study. 

VI. CONCLUSIONS 

We present three-dimensional realizations of the KBC 
model for the D3Q27 lattice. We review the details of 
the entropic stabilization and describe eight variations of 
KBC. Stability and accuracy is studied in detail for ho¬ 
mogeneous isotropic turbulence. A detailed comparison 
with LBGK is carried out at various grid resolutions. Sec¬ 
ond order rate of convergence is numerically confirmed in 
the vast majority of the statistical quantities of interest. 
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It must be stressed that the entropic KBC models are 
stable (in contrast to LBGK and RLB) for all the con¬ 
sidered cases here; despite of under-resolution and high 
Reynolds numbers. 

The KBC models were shown to capture the expected 
scaling law for energy spectra in the case of high Reynolds 
numbers. Low order statistics such as averages of kinetic 
energy, enstrophy and rate of dissipation as well as the 
spectral densities for energy and rate of dissipation agree 
well with resolved simulation despite under-resolution. 

In general, we show that by keeping the kine¬ 
matic (shear) viscosity coefficient constant the presented 
method is extremely stable and produces accurate re¬ 
sults in presence of under-resolution. These findings and 
the parameter-free and explicit nature of KBC as well as 
the lack of explicit turbulence models renders the scheme 


a very promising candidate for applications in both re¬ 
search and engineering contexts were high Reynolds num¬ 
bers and computational cost are of importance. 

It has been demonstrated in [T] that also for three di¬ 
mensional flows in presence of complex walls low-order 
statistics can be captured well using KBC. In a further 
publication we will address the issue of boundary condi¬ 
tions for KBC models in both two and three dimensions. 
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Appendix: Central Moments 

The mapping between natural and central moments is linear in the non-conserved moments and given by the 
following relations: 




^xz — ^xz '^x'Uz 


ffys — ^yz Uyllz 


Nxz = -ul+ul 


^yz — ^yz '^y '^z 


T — T — {u^ Uy u^) 


Qxyz — Qxyz '^x^^yz '^y^^xz '^z^^xy Ux^y'^z 


Qxyy — Qxyy 3 (^'^y^xy H“ Ux{^Uy + 2 Nyz ^xz 


Qxzz = Qxzz - 5 {^ujlxz + Uxiiul - Na;z - Nyz + T 


Qxxy — Qxxy 3 (^'^x^xy “b ^yi.^'^x ~b ‘^^xz ^yz “b 


(A.IO) 


Qyzz = Qyzz “ 5 {^ujly:, + Uy{3ul - - Ny^ + T 


(A.ll) 


Qxxz — Qxxz ~ h { (iUx^xz + Uz{3u^ + 2N^z ~ Nyz + T 


(A.12) 


Qyyz — Qyyz ~ | {^UyYlyz + Uz{3Uy + 2Nyz — Nxz + 


(A.13) 


Mq22 — Mq 22 — \^Uyiil + 4uyUzTlyz + (Uy + u1)T/3 — (Uy + u‘l)Nxz/3 + (—Uy + 2u1)Nyz/3 

“b 2UyQyzz “b 2UzQyyz'\ 


(A.14) 


M 202 — M 202 ~ + ‘^UxUzTlxz + + u1)T/3 + (—+ 2u‘l)Nxz/3 — (u^ + u1)Nyz/3 

-\-2uxQ xzz + 2uzQ xxz\ 


(A.15) 


M 220 — ^220 ~ + '^UxUyllxy + {u^ + Uy)Tj3 + {—u^ H- 2uy)Nxzl^ + — Uy)Nyz/3 

-\-2u 

xQxyy H” 2,y^yQxxy^ 


(A.16) 
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-^^211 — ^211 ^^j.1lyUz ^UxUz^xy ^Ux'^y^xz ^x^yz “t“ UyUzT/“^ “h ‘2'liyUzNxz/'^ '^y'^z^yz 

~\~^^xQxyz “1“ ^zQxxy ^yQxxz i 


(A.17) 


Mi2l — Mi2l ^'Ua^Ziy'U^ “I” ^UyUz^xy ~\~ UyTlxz H“ ^UxUy^yz H“ UxUzTj^ UxUzNxz H“ ^UxUz^yz/^ 

-\-^UyQxyz H“ '^zQxyy H“ '^xQyy::^ 


(A.18) 


Afll2 — A/ii 2 “1“ -h ‘lUxUz^^yz “t“ '^a:'t^y7^/3 '^x'^y^xz'^x'^y^yz 


~\~‘^UzQxyz H“ UyQ xzz “1“ Q 


y 22 


(A.19) 


Mi22 = Mi22 - (u^ulul + 2uyulll^y + 2ulu^Il^;, + Au^UyUjiy:, + {u^u^ + u^ul)T/i 


+ ( Ux'lly '^x‘^z')^XZ l‘^ ( '^X^y “t“ ‘^UxUz)^y^ 1^ 

~\~^UyUzQxyz “1“ '^zQxyy '^yQxzz H" ‘^'^x'^yQyzz ^'^x'^zQ 

-\-2UyMii2 + 2UzMi2l + Ux^022 


yyz 


(A.20) 


^^212 — A7212 “h 4zia;Wy'U2na;2 A 2w^'U^ny2: 4“ (^x^y '^y^z')'^!^ 


-\-{-ulUy + 2UyUl)Nxzl^ + (-'^X% - '^y'^l)^yzl^ 


~\~^'llxUzQxyz ‘^^x'^yQxzz '^zQxxy 

-\-2UzM2ll + 2 Wa;Afii 2 + UyM202 


■ U^QyzZ 


2'lly'lizQx 


(A.21) 


A /221 — Af221 ^Uj.UyUz ~\~ '^Ux^y'^z^xy “t“ 2'Ua;'UyIIa^2: “t“ “t“ {,’^x^z 4“ 

+(-^x^^ + 2uluz)Nxzl'^ + - u^ynz)Nyzl‘^ 

~\-4:UxUyQxyz “t“ 2UxUzQxyy “t“ 2UyUzQxxy “t“ UyQx 

+2UyM2ll + 2 UxMi2l + UzM220 


■ '^xQyyz. 


(A.22) 


M 222 = Ar 222 - (ululul + AuxUyulUxy + AuxuluzUxz + 4u^? 


+(—u^u 


+ 4u^%Uznj,2 + (u^My + + UyUl)T/i 

■ V + 2ulul)Nxz/3 + i-ulul + 2ulul - ulul)Nyz/3 

~\~SUxUyUzQxyz A 2UxUzQxyy A 2UxUyQxzz A 2Uyll^Qxxy A 2Vj^UyQyzz A 2UyUzQxxz A 2zi^ 

+U^Mq22 + U^M202 + u\M220 + 4UyM2M2ii + AUxUzMi2i + AUxUyMii2 + 2UxMi22 + 2UyM2l2 + 2UzM'^ 


luzQ 


yyz 


221 


(A.23) 


Substituting Eqs. (A.1)-(A.23) in Eq. (19) one arrives at the moment representation in the central basis: 
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f(0,0,0) = pi- (ul - 1 ) (ul - 1 ) (ul - 1 ) 

(^2 l) Ali^Uz i^y ^xz ^UyU^; (^ 3 ; l) z 

+ (“x {’’^y + “z ~ 2 ) + (1 — 2Uy) + Uy) -/Va;z/3 + (Uy [u^. + — 2) + (l — 2^^) + M^) Ny^ 


■ (ulul 4 

2 2 1 2 2 

- M^M^ + MyM^ — 

Hui 

+ ul+ul) 

+ 3) r/3 

00 

Qxyz 





2ux (ul 

1) Qxyy 

2ux 

i<- 

1) Qxzz 


2Uy iul 

1) Qxxy 

2Uy 

(ul- 

1) Qyzz 


2uz [ul 

1) Qxxz 

2uz 

(ul- 

- 1) Qyyz 


-(l-< 

) -^022 + (1 

-< 

;) M202 + (1 — ^2 

) M220 


-AUyUzM211 - ^UxUzMi 21 - AUxUyMii2 
— 2UxMi 22 — 2UyM2l2 — 2 UzM221 — ^ 222 ^ 


f{a,0,0) = Ip (SUa; (ul - l) (m^ - l) (ct + Ua;) 

+6Uy (ul - l) (cr + 2Ma;)na;y + [Uy - l) (cT + 2Ux)^xz + 12Ma;%Uz(o- + Ux)^yz 

- [uux (Wy + ul-2) +ul [ul + - 2) - 2 (uy - 1) [ul - l)) 

+ [ul i2Ux{(T + Ux) - ul) - Ux (ul + 1) (cr + Ux) + My + - l) Ny^ 

+ (ul (uxicr + Ux) + ul) + Ux {ul - 2) (ct + Ux) - My - m^ + l) T 
-t-12MyM2(cr -t- 2Ux^Qxyz 

+3 (ul — l) (cr + 2Ux)Qxyy + 3 (m^ — l) (ct + 2Ux)Qxzz 
+6Uy (ul - l) Qxxy + 6UxUy{a + Ux)Qyzz 
+6 m 2 iUy Qxxz “t” ^UxUzic Ux)Qyyz 

+3Ma;(cr + Ux)Mq22 + 3 (m^ — l) M202 + 3 (m^ — l) M220 

+ 12MyM2-M211 + 6Uz{o’ + 2Ux)Mi2i + 6My(cr + 2Ma;)-Mii2 

+3 (cr + 2Ma;) M 122 + 6MyiW212 + GUzM 221 + 3i\cf222^ 

f{0,\,0) = Ip (3My (ul - 1 ) (ul - 1 ) (A + My) 

“t“6Ma; (m^ (A “t“ 2My)na:y “t“ ^2UxUyUz{,^ Uy^TYxz “t“ Sm^ (m^ (A “t“ 2Uy')Tl.yz 

+ { — '^X ('^y(A + Uy) + ul — 1 ) + Uy (2M^ — 1 ) (a + Uy) + M^ — l) Nxz 
+ i — ul ^My(A + My) — 2ul “1” 2^ — My ^M^ — 2^ (A + My) — 2 m^ + 2^ ^yz 

+ iux (My(A + My) -t- M^; - ^y i^Z - 2) (A + My) - U^ + 1 ^ ^ 

-j-12Ma;M2 (A “1- 2Uy)Qxyz 

“l“6Ma; (m^ Qxyy “t“ ^UxUy{\ -\- Uy)Qxzz 

+ 3 (m^ — 1 ) (a + 2Uy)Qxxy + 3 (m^ — 1 ) (a + 2Uy)QyzZ 

~\~QUyUz{X Uy)Qxxz “t” ^Mg ^M^ Qyyz 

+3 (m^ — 1 ) Mq22 + 3My(A + Uy)M202 + 3 (m^ — 1 ) M 220 
+6 m 2(A + 2My)Af211 + 12Ma:M2Afi21 + 6Ma;(A + 2My)iWil2 
+6Ma;-/Vfi22 + 3 (A + 2Uy) M 212 + 6MZ-M22I + 3A/222^ 


(A.24) 


(A.25) 


(A.26) 
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/( 0 . 0 . 5 ) = Ip (3m^ {ul - 1) {ul - 1) ((5 + U^) 

-\-\‘2'Ujx’^yU 4“ ^z}^xy 4“ (^y 4“ ‘^‘^z')^xz 4“ 6?Xy ^^2; 4“ 2'U2:)ny2; 

4" {~'^x (^z(^ 4- Uz) 4“ Uy — l) + Uy{2uz(6 + Uz) 4- 1) — Mz('^ 4- Uz) ~ l) 

+ (ul (2Uz(S + Uz) -u^ + l) - ul(UziS + Uz) - 1) - u^{5 + Uz) - l) Nyz 

4" (^z(^ 4“ Uz) 4“ My — l) 4- ul(uz(5 4- Uz) ~ 1) ~ 2uz(S + Uz) 4- l) T 

-l-12Ma^My ((5 “t" 2Uz')Qxyz 

-\-QUxUz(S + Uz^Qxyy 4“ ()Ux (Uy Qxzz 
+ 6 MyMz(^ 4“ Uz)Qxxy 4“ 6 My (Ux l) Qyzz 
+ 3 (Uy — l) (5 + 2Uz)QxXZ 4" 3 (m^ — l) (<5 + 2Uz)Qyyz 
4"3 (ul — l) Mq 22 + 3 (uy — l) M 202 4- 3uz(S + Uz)M 22 o 
+6uy(S + 2 uz)M2ii 4- 6ux(6 + 2uz)Mi2i 4- l2uxUyMii2 
-\- 6 uxMi 22 4" 6 uyikf 2 i 2 4- 3 ((5 + 2uz) M 221 4- 3 Af 222 ^ 

/(^,A.O) = \p (-UxUy (ul - 1) (A + My) (cr + Ux) 

(Uz 4“ 2My)(fJ 2Ma;)na;y 2MyMz(A “t" Uy)((J “t" 2M2;)n2;Z UxUz(^ 4“ ‘^Uy) (<7 “t" M^;) Ily z 

+ 1 (aUx (Uy(X + Uy) +ul-l) +ul (Uy(X + Uy) + M^ - l) - 2Uy (ul - l) (A + My)) lVa;z 

+ ^ (cUx (Uy(X + My) - 2M^ + 2) + ul (Uy(X + My) - 2M^ + 2) + Uy (ul - l) (A + My)) Nyz 

+ 1 ( — aUx (Uy(X + Uy) + M^ — l) — M^ (My(A + My) + M^ ” 1) ” My (m^ — l) (A + My)) T 

2Mz(A-1- 2Uy)(o' 2Ux)Qxyz 

(^z ^) 4“ 2Ux)Qxyy Uy(X + Uy)((7 + 2Ux)Qxzz 

— (ul - 1) (A + 2My)Q 

a:a:y (A + 2'Uy)((7 + Ux)Qyzz 

-2UyUz{X + Wy)Q a:a:2: 2UxUz(( 7 -\- Ux)Qyyz 
— Ux((J + Ux)Mo 22 — Uy(X + My)M202 4" (l — M^) M 22 O 
-2mz(A + 2My)M2ii - 2Mz(cr 4“ 2ux)Mi2i - (A + 2uy)(a + 2ux)Mii2 

— ((7 + 2Ux)Mi22 — (a + 2Uy)M2l2 ~ 2Mz-M 221 ~ Af222^ 

/(<t. 0,5) = jP (-UxUz (ul - 1 ) (5 + Mz)(cr + M^) 

2UyUz(5 4“ Uz)((J 4“ 2Ma^)na;y (^y ^) (*^ 4“ 2Uz)(o' 2Ma^)na;z 2UxUy(6 -j- 2Uz)(o' -t- Ma;)IIyz 

+ 1 (ctM^; (Uz(6 + Uz) + My - 1) + M^ (Mz((5 + Mz) + My - 1) - 2 (m^ - l) Mz((5 + Mz)) Nxz 
4-| (crUx ( — 2Uz(6 + Mz) 4- My — 1 ) + M^ ( — 2Uz(S + Uz) 4- My — 1 ) + (My — 1 ) Uz(6 + Uz)) Nyz 
+ 1 (-CrUx (Uz(S + Uz) + My - 1 ) - M^ (Uz(S + Uz) + My - l) - (Uy - 1 ) Uz(S + Uz)) f 
2Uy(S 2Mz)(c 7 “t" 2Ux)Qxyz 

Uz(^ 4“ Uz)((y 4“ 2Ux)Qxyy (Uy 1 ) (fT + 2Ux)Qxzz 

2UyUz(S Uz)Qxxy 2UxUy((7 -j- Ux)Qyzz 
(Uy 4 ) 4“ 2Uz)Qxxz Ux(S 4“ 2Mz)(<7 + Ux)Qyyz 

-Ux((7 + Ux)Mo22 4- (1 - My) M 202 -Uz(6 + Uz)M220 

— 2Uy(S + 2Mz)Af211 — (|5 4“ 2Uz)(<J 4“ 2Ux)Ml21 ~ 2Uy((T + 2Ma;)Afii2 

— (cr + 2 Ux)Mi 22 — 2 UyM 2 l 2 — ((5 + 2 Uz)M 221 ~ -Af 222 ^ 


(A.27) 


(A.28) 


(A.29) 
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f{0,\,s) = Ip {-UyU^ {ul - 1) ((5 + u^){\ + Uy) 

‘^UxUyi^^ ~\~ 2'Ujg)(A -t- Uy^^xz {^x ^) (*^ ^~ 2zi_g)(A “t- 2Wy)ny_g 

+ ^ (AUy (— 2Uz{S + Uz) + (— 2Mz( 5 + Mz) + — l) '^z{^ + ^z)) 

+ 1 (AUy (uzi6 + Uz) +ul-l) +ul (uz(S + Uz) + U^ - l) - 2 (u^ - l) Uz(S + Uz)) Nyz 
+ i {-Xuy {uz{6 + Uz) + ul- l) - Uy (uzi6 + Uz) + - l) - {ul - l) Uzi6 + Uz)) f 

2Ux{^ 2Uz){X -j- 2Uy)Qxyz 

~ ~ (A.30 

2UxUz{S Uz)Qxyy 2iUxUy{\ A Uy)Qxzz 

Wz(^ A Uz){X 2Uy)QxXy {'^X A 2Uy)Qyzz 

2Uz){X Uy)QxxZ {'^x A 2Uz)Qyyz 

+ (l ~ ’^x) Mo 22 — Uy{\ + Uy)M202 — Uz{5 + Uz)M220 

— {6 + 2uz){X + 2uy)M2ii — 2,Ux{S + 2uz)M\2i — 2ux{X + 2uy)Mii2 

— 2UxMi22 — (a + 2Uy)M2l2 — ((5 + 2Uz)M221 — M222^ 

f{cr,\,5) = Ip [UxUyUziS + Wz)(A + Uy){a + Ux) 

-\-Uz(^S 'Uz)(A 2Uy))cr 2Ux)^xy A A ‘2Uz){X A Uy))(T A 2Ux)^xz “t“ ^xi^ “t“ 2'Uz)('^ “t“ A ^x)^yz 

a| {-Suz{ux{(J a Ux) - 2uy{X A Uy)) - ul{ux{(T A Ux) - 2uy{X A Uy)) - UxUy{X A Uy){a A Ux)) Nxz 
A^ {^^z i‘^^x ~t“ ^x) ^y)X A ’^y)) A U^{2Ux{jy A Ux) ^y(.X A Uy)) UxUy)X A ^y)(^ A Ux)) XJ^yz 

A^ (Suz(uy(X A Uy) A Uxia A Ux)) A ul{uy{X A Uy) A Ux{<7 A Ux)) A UxUy{X A Uy){a A m^)) T 

A((5 a 2uz){X a 2uy){(j A 2ux)Qxyz 

~\~Uz{5 A Mz)(o' A 2Ux)Qxyy A Uy{X A Uy){(T A 2,Ux)Qxzz 
A'Uz(^ A Uz)(^X A 2uy)Qxxy “t“ ^x{.X A ‘2uy))o' A Ux)Qyzz 

-\-Uy{6 A 2uz){X A Uy)Q XXZ H“ (J H- 2'Uj2)((J H- Ux^Qyyz 
-\-Ux{<J + '^a:)-^022 + Wy(A + Wy)M202 + + '^^)A^220 

+((5 + 2'U2:)(A + 2uy)M2ii + (^5 + 2uz){o' H- 2'Ux)-A^i2i + (A H- 2uy){(j H- 2ux)^ii2 

+ ((7 + 2'Ux)Afi22 H“ (A + 2Uy)M2l2 H“ (<5 + 2Uz)M221 H“ ^222^ 


(A.31) 



